A Neural-Network-Optimized Hydrogen Peroxide Pairwise Additive Model for Classical Simulations

We have developed an all-atom pairwise additive model for hydrogen peroxide using an optimization procedure based on artificial neural networks (ANNs). The model is based on experimental molecular geometry and includes a dihedral potential that hinders the cis-type configuration and allows for crossing the trans one, defined between the planes that have the two oxygen atoms and each hydrogen. The model’s parametrization is achieved by training simple ANNs to minimize a target function that measures the differences between various thermodynamic and transport properties and the corresponding experimental values. Finally, we evaluated a range of properties for the optimized model and its mixtures with SPC/E water, including bulk-liquid properties (density, thermal expansion coefficient, adiabatic compressibility, etc.) and properties of systems at equilibrium (vapor and liquid density, vapor pressure and composition, surface tension, etc.). Overall, we obtained good agreement with experimental data.


■ INTRODUCTION
Hydrogen peroxide 1 is the simplest of the peroxides and the smallest molecule that exhibits enantiomeric isomerism, making it a chiral compound. 2Its pure liquid state is denser than water but only slightly more viscous.Like water, it forms hydrogen bonds, which give it a wide temperature window for the liquid state.In fact, the difference between its boiling point, 423.3 K, and melting point, 272.72 K, 3 exceeds that of water, and its vapor pressure is indeed lower, 5 mmHg at 303 K. 4 Given its chemical similarity to water, it is not surprising that the liquids are miscible for all concentrations and temperatures.It is worth mentioning that it produces a eutectic mixture with water showing a strong dependent melting point with composition. 4here are several uses for hydrogen peroxide in daily life and in industry.When diluted in water, it is used as a bleaching agent, including in hair bleach, and as an antiseptic. 5It is commonly used at home to disinfect sores and is part of some products to whiten teeth. 6In industry, it is often used at higher concentrations, and most notably, it can be used as a rocket propellant. 7These applications are based on the fact that peroxide is unstable and an oxidizing agent.
Its molecular structure has two oxygen atoms joined by a single bond, each bonded to its corresponding hydrogen atom (see Figure 1).Thus, its chemical formula is H 2 O 2 .The molecule is not planar and at its configuration of minimum energy shows a certain dihedral angle, defined by the planes containing the oxygen atoms and each hydrogen, close to 110°.This configuration is compatible with two enantiomers, structures that are mirror images of each other, which can be transformed one into the other by crossing a small potential barrier associated with the trans (planar) configuration. 8,9In principle, the enantiomers are also connected through the cis (also planar) configuration, but the height of the energetic barrier is much larger through this path. 8,9lthough interesting and small, the hydrogen peroxide molecule has not been extensively studied by simulation techniques.−18 In the latter lists, the classical models given by Chung 16 were employed to compute residence times in crowded protein environments and diffusion on the surface of a protein.−21 From the possible models, we would like to highlight the one recently given by Orabi and English. 18The model is based on the CHARMM force field, 22,23 where the geometrical parameters are taken from an ab initio geometry optimization of the H 2 O 2 single molecule, the partial charges are adjusted to match the dipole moment and energy as a function of the dihedral angle (also from ab initio computations), the hydrogen Lennard-Jones (LJ) parameters correspond to the TIP3P model 24 (CHARMM force field 22,23 ), and the oxygen LJ parameters are fitted to reproduce the density and heat of vaporization at 0 °C.That is, the authors of this model followed a well-accepted and popular path to yield their H 2 O 2 model (there are other ways such as fitting the relative dielectric constant and the temperature of maximum density, 25 for instance).However, to our knowledge, there are still no simulations providing thermodynamic data such as boiling point, vapor pressure, surface tension, critical point, and viscosity, among other important properties for pure peroxide and water solutions.
In this study, we propose an alternative model for H 2 O 2 based on the OPLS-aa force field 26−28 to address the lack of an intensively tested H 2 O 2 model.Note that, unlike the CHARMM force field, the OPLS-aa does not include LJ parameters for the H sites, simplifying the model.We parametrize the pairwise additive model using artificial neural networks (ANNs) trained on corresponding sets of molecular dynamics results.To do so, we define ranges for the model's geometrical structure, including bond distances and angles, LJ parameters of the oxygen atoms, partial charges, and the energy contribution associated with the H 2 O 2 dihedral angle.We use the molecular geometry of the solid and vapor phases as a basis, but, in contrast to the method used by Orabi and English, we sweep all possible parameter values into ranges to minimize differences between the ANNs' predictions and experimental data.Our approach aims to optimize the parameters by considering several properties of pure H 2 O 2 and water mixtures simultaneously, targeting 22 different values.−41 We present several liquid and vapor properties of the obtained pairwise additive model as a function of the weight percentage of the water− peroxide mixture and temperature.

■ MODEL
The geometrical parameters necessary to define the model are the oxygen−oxygen distance, d OO , the oxygen−hydrogen distance, d OH , the angle formed between the lines joining the oxygen−hydrogen and oxygen−oxygen sites, a HOO , and the dihedral angle given by the position of all atoms, a d .These parameters and the corresponding force constants define the so-called bonded interactions.We use the OPLS-aa force field 26−28 potential functions, including a Ryckaert-Bellemans (RB) dihedral contribution.Thus, we employ harmonic potentials for bonds and angles and the dihedral function  2).Note that the large and tiny barriers correspond to the cis and trans configurations, respectively.c 2 controls the height of the large barrier, which translates into the position of the a d peaks.Thus, the selected range of c 2 is such that the a d peaks are located close to the experimental values for the solid and vapor phases. 8,42The so-called nonbonded interactions are given by the LJ potential acting only on the oxygen sites, with parameters ϵ OO and σ OO , and the Coulomb potential acting on every site.We are excluding all internal LJ contributions to the model.The values of the partial charges located at hydrogen sites, q H , and oxygen sites, −q H , define the intra and intermolecule Coulomb interactions.
The geometrical parameters corresponding to the solid 42 and vapor 8 phases are 0.1453 nm, 0.0988 nm, 101.9°, 90.2°, and 0.1474 nm, 0.0970 nm, 94.8°, 115.5°, corresponding to d OO , d OH , a HOO , and a d , respectively.Thus, we are setting the ranges for the geometrical parameters in 0.145−0.148nm for d OO , 0.094−0.099nm for d OH , and 94−102°for a HOO .This way, we cover the solid−vapor experimental range while still allowing for some flexibility in the parameters.For the case of a d , as mentioned in the previous paragraph, the selected range of c 2 combined with the other parameters yields a d values inside the experimental range.In addition, the q H range found in the literature is very large, 11,19−21 and we need to narrow it

Journal of Chemical Theory and Computation
down.We have noticed that the Atomic Charge Calculator II, 43 using empirical methods along with parameters from literature, leads to values of q H between 2 and 10% larger for H 2 O 2 than for H 2 O when employing the same method.Also, Martins-Costa and Ruiz-Loṕez yields similar results. 11In addition, three-point water models such as SPC, 44 SPC/E, 45 and TIP3P 24 have q H values in the 0.41−0.4238e range.Hence, we find it reasonable to set the q H range in 0.43−0.47e.We have also noticed that this range yields relatively good results when combined with the other parameter ranges.Finally, we observe that larger partial charges produce very large viscosities.We set the σ OO range in 0.296−0.300nm to get densities close to the experimental ones, since the impact of σ OO is strong on this property.The ϵ OO range in 0.7−0.9−28  The model itp file, defining the optimized parameters and the force constants, is included in the first section of the Supporting Information (SI).The optimized parameters are d OO = 0.1463 nm, d OH = 0.0979 nm, a HOO = 95.79°,q H = 0.4323 e, σ OO = 0.29989 nm, ϵ OO = 0.8912 kJ mol −1 , and c 2 = 26.91 kJ mol −1 .In the following sections, we describe the route we followed to reach these optimized values.
Simulations Details.We employ the gromacs 2022.4 molecular dynamics package 47−50 to produce and analyze trajectories of H 2 O 2 −H 2 O mixtures.In all cases, we use the SPC/E rigid water model. 45This model has a good performance considering it has only three sites. 51,52The weight percentage of H 2 O 2 , wt %, in the solutions is 0, 10, 20, 30, 40, 60, 80, and 100.The number of H 2 O 2 and H 2 O molecules for each box is given in such a way as to yield similar cell sizes at room conditions and a minimum of 6000 oxygen atoms.Hence, we use 6000 H 2 O molecules for 0 wt % and 4350 molecules for 100 wt %.Note that the number of oxygen atoms increases with wt % due to the higher density of H 2 O 2 .Before conducting the formal runs described below, we built the system and performed the necessary tasks to achieve thermalized initial configurations.
To access densities, ρ, the relative dielectric constant, ϵ, diffusion coefficients, D, enthalpy of vapor formation, ΔH vap , coefficient of thermal expansion, α P , isothermal and adiabatic compressibilities, κ T and κ S , number of hydrogen bonds, n h and structural properties, we produce 100 ns trajectories with 0.5 fs time steps for bulk-liquid systems.For properties derived from fluctuations, 53 we carry out a block-averaging error analysis and discard those blocks far (more than three standard deviations apart) from the average.We discarded only three blocks in all runs for the determination of ϵ, all placed at their beginning (the thermalization length was probably not enough in these cases).The trajectories are written with a frequency of one each of 1000 steps.For this purpose, we set the v-rescale algorithm with a 0.1 ps coupling time for temperature coupling and the c-rescale isotropic algorithm with a 1.0 ps coupling time for pressure coupling.We set a 1.2 nm cutoff length for nonbonded interactions, considered dispersion corrections for energy and pressure, and used a Fourier spacing of 0.12 nm to handle the long-range Coulomb interactions with the PME method.The properties are obtained by using the energy, dipoles, msd, hbond, and rdf gromacs tools with the proper flags.These relatively large runs are needed to avoid large uncertainties for the relative dielectric constant.
Shear viscosities, μ, are accessed through independent 50 ns trajectories, setting the amplitude of the acceleration profile to 0.05 nm ps −2 (nonequilibrium dynamics) and turning off the pressure coupling. 54For these runs, we have also turned on the constraints to h-angles, meaning that d OH and a HOO get fixed but not the dihedral angle nor the d OO distance.By doing so, we get results much closer to the experimental values reported by Phibbs and Giguere 55 (around 30% smaller).For all the explored values of the parameters, the model most frequently produces viscosity values above the experimental ones 55 for T = 293 K.
Finally, we access the surface tension, λ, vapor pressure, P v , vapor composition, wt % v , boiling temperature, T b , critical point temperature, T c , pressure, P c , and density, ρ c , by setting 1 × 1 × 3 aspect ratio prismatic simulation cells containing a liquid slab in equilibrium with its corresponding vapor phase. 56,57In this case, we sampled from the NVT ensemble, turned off the energy and pressure dispersion corrections (the system is not homogeneous), increased the cutoff distance to 1.6 nm (to partially compensate for turning off the dispersion corrections), and set a number of Fourier wavevectors along the largest side of the prismatic cell 3 times larger than those of the other two sides.In addition to these details, the system size must be sufficiently large to avoid spurious effects 58 (the minimum of 6000 oxygen atoms is well-above the limit where these undesired effects appear).
The above-given paragraphs describe the simulation details to produce the formal runs for the final results of the model.However, since we need a large number of determinations to carry out the parametrization, we performed much shorter runs (20 ns with a time step of 1.0 fs) with smaller system sizes for this purpose.To study bulk-liquid systems, we consider 1424 peroxide and 1152 water molecules to yield a 70 wt % peroxide−water mixture and 2000 peroxide molecules for pureperoxide.To access the liquid−vapor properties, we set 2136 peroxide and 1728 water molecules and runs of only 10 ns.From these simulations, we get ρ, ϵ, μ, P v , λ, and wt % v , for pure and wt % = 70 mixtures at different temperatures and as a function of the varying parameters to train and validate the ANNs (see Table S1 and S2 of the SI).
Model Parametrization.We produced 1000 determinations (Table S1 of the SI) of the target properties to train the ANNs.We normalized the input and output vectors (composed of 7 and 22 items each) to [0,1] and employed the backpropagation algorithm.The target properties correspond to pure H 2 O 2 and 70 wt % of water−peroxide mixtures.We considered the evaluation of ρ, ϵ, μ, P v , λ, and wt % v .The last property applies only with the 70 wt % mixture.The corresponding experimental data and their sources are given in Table S1 at the SI.
Once trained, the predictions of the ANNs are compared with another set of 300 determinations of the target properties to evaluate their behavior.These 300 determinations (Table S2 of the SI) are obtained following the same procedure as for Table S1.We tried sequential ANNs having from 1 to 3 layers, with 2, 4, 8, and 16 neurons per layer, N L .We observed that

Journal of Chemical Theory and Computation
further increasing N L worsens the results.Also, we tried ReLu and Sigmoid functions following the linear steps to introduce nonlinearity.We use the PyTorch framework 59 to implement the ANNs.Table S3 of the SI shows the deviations of their predictions from Table S2.Arguably, the best results are observed for ReLu/1/8, ReLu/1/16, ReLu/2/4, ReLu/2/8, Sigm/1/8, and Sigm/2/4 ANNs, where we follow the notation nonlinear-function/number of layers/N L .Note that a relatively small number of neurons can capture the general behavior of the data, something not easily observed by the naked eye.We believe that the relatively large dispersion of the simulation data, due to the lack of sufficient statistics (short runs), is smoothed out by the ANN training process.Note that the number of parameters, N P , increases with N L and the number of layers.Thus, a large number of neurons tend to capture not only the genuine trend but also the noise of the training set, which negatively impacts its predictions.
As mentioned, the seven input parameters, p, sweep the following data ranges: 25.5 kJ mol −1 ≤ c 2 ≤ 36.0 kJ mol −1 , 0.43 e ≤ q H ≤ 0.47 e, 0.296 nm ≤ σ OO ≤ 0.300 nm, 0.7 kJ mol −1 ≤ ϵ OO ≤ 0.9 kJ mol −1 , 0.145 nm ≤ d OO ≤ 0.148 nm, 0.094 nm ≤ d OH ≤ 0.099 nm, and 94°≤ a HOO ≤ 102°.We hoped that, with these ranges, a clear local minimum of a target function could be found far from the borders of the ranges.We define the following target function where x i is the value obtained for the property at given conditions, x ei is its corresponding experimental value, and i runs on all (22) properties and conditions.In general, wc i is 1 for the pure substance and 0.25 for the mixture with the SPC/ E model, wp i is 1 for ρ, ϵ, P v , and λ and 0.5 for μ.The rationale behind our wc i selection is to incorporate the water−peroxide interaction in the training process while avoiding a large dependence of the specific choice of a water model.In addition, we have decreased some weights where we have observed strong fluctuations from the simulation results (the second row of Table S5 at the SI shows the weights).The weights are somewhat arbitrary, but the idea is to obtain a holistic model capable of capturing several properties at the same time.The same procedure with other weights could lead to a model able to better approach some particular properties of the list.For all trained ANNs, we generate 1 × 10 6 random p vectors to make predictions that are evaluated through eq 1.From this procedure, it turns clear the existence of several F(p) local minima for all ANNs.Next, we select the p reaching the F(p) minimum and locally search for even lower values by slightly changing p. Table S4 of the SI shows the final outcomes.Note that some but not all p from different ANNs are similar.Again, this occurs since there are several combinations of parameters leading to similar F(p) values.Note also that some p are very close to the borders of the defined ranges.This may suggest that a better fit to the experimental data could be obtained outside the boundaries we have imposed.However, it should be noted that for all of these cases, the predicted a HOO values were either 94 or 102 degrees, which are beyond the experimental determinations for the solid and vapor phases.
Finally, starting from the values given in Table S4 that are not at the border of the ranges (we avoid getting outside the ranges), p start , we perform a local search by carrying out more MD simulations with p try = p start + ξ ⊙ Δp, where ξ is a vector of [−0.5,0.5]homogeneously distributed random numbers, ⊙ is the Hadamard product, and Δp has equal nonzero components in normalized space and a small magnitude.We present the best six p try in Table S5.It should be mentioned that all these p try produce low F(p try ) with all trained ANNs.From this list, we have selected the p try from the first row, i.e., from the model with the smallest error.However, due to fluctuations in the results and given the similar values of the errors, any of the six combinations should behave similarly.Note also that rows 1 and 6 and 2−5 are akin.

■ RESULTS
Structural Properties of the H 2 O 2 Bulk Liquid.For pure H 2 O 2 , T = 293 K, P = 1 bar, and by following the procedure described in the second paragraph of the Simulations Details section, we build Figure 2. The figure depicts the most important structural properties obtained for the bulk-liquid phase.These are the radial distribution functions for all sites (obtained by using the rdf tool), the probability density functions corresponding to the distribution of angles a HOO and a d (angle tool), and the dipole moment distribution (dipoles tool).We have also added the Ryckaert-Bellemans (RB) potential to show its effect on the a d distribution.
The O−O g(r) shows a broad main peak around 0.30 nm, followed by a valley and a second-shell shallow maximum.Indeed, the peak shape, its height, and the complete g(r) are similar to those shown by Orabi and English, 18 although the temperature is T = 273 K in their case.In contrast, our g(r) strongly differs from the one shown by Yu et al. 60 The broad main peak appears to be the result of two overlapping peaks, indicating a parallel configuration of oxygen atoms where the two pairs are positioned closely to each other.Therefore, the cross-neighbors contribute to a peak that appears slightly further away.This finding aligns with the two most stable conformers of two-molecule aggregates observed in reference 18, as shown in Figure 2 of that reference for snapshots of several aggregates in the vapor phase of the pure substance.The H−H g(r) is also similar to the one shown by Orabi and English, 18 having two close peaks at 0.23 and 0.27 nm.The former is due to the intramolecular contribution, whereas the second is due to the H−H distance involving a hydrogen bond. 18Again, the general shape, locations, and heights of the peaks are much closer to those given by Orabi and English than from the work of Yu et al. 60 However, the principal peak of our O−H g(r), the one corresponding to the hydrogen bonds, has a height close to the value given by Yu et al. (around 4.5), contrasting with the smaller value from Orabi and English (around 3.5).The location of this sharp peak is close to 0.188.Finally, the general shape of the O−H g(r) resembles the one presented by by Orabi and English. 18he a HOO probability distribution function, pdf, is given in the inset of Figure 2a.The peak of the distribution is placed around 96°, in good agreement with the minima set for this angle.In addition, it is closer to the value of the vapor phase than the one of the solid, close to the values set by Caputo et al. 13 and to the experimental values given by Redington et al. 61 Figure 2b depicts the dihedral a d pdf and the RB contribution to the energy, both even functions peaking at ±104 and 0°, respectively.The peak of a d is close to the one found by Moin et al. 12 Note that there are no cis configurations, given the large RB peak, but trans configurations are allowed.For our model, the trans state turns a 100% pure right enantiomer liquid in a 50−50% mixture in a few tens of picoseconds.Finally, the inset of Figure 2b shows the dipole pdf of H 2 O 2 molecules.The distribution peaks at 2.69 D and is asymmetric, showing a tail toward zero that correspond to the trans configuration.Note that μ monotonously increases for decreasing |a d |.
Water−H 2 O 2 Bulk-Liquid Mixtures.Figure 3 shows several properties of the bulk-liquid mixtures of H 2 O 2 and H 2 O as a function of the weight percentage of H 2 O 2 , wt %, and T.These are the density, ρ, the thermal expansion coefficient, α P , the adiabatic compressibility, κ S , the relative dielectric constant, ϵ, the enthalpy of vaporization, ΔH v , the diffusion coefficient of H 2 O 2 and H 2 O, D, the viscosity, μ, and the number of hydrogen bonds per oxygen atom, n h .
As expected, ρ increases with wt % and decreases with T. Furthermore, the ρ increasing rate with wt % also augments with wt %.The inset of panel a of Figure 3 compares the simulation results with experimental data from Easton et al. 62 for T = 273 K (circles and black line) and T = 313 K (triangles-left and cyan line).The agreement is very good.In addition, the obtained results for the SPC/E water model agree well with those reported in previous works. 52The ρ dependency with temperature allows us to compute by fitting a cubic polynomial to ρ(T) from where we obtain its partial derivative appearing at the right-hand side of eq 2. α P can also be computed from fluctuations 53 by using the energy tool, which leads to a general good agreement with the fitting procedure.However, we have noticed that the fitting procedure yields smoother trends, the ones reported in Figure 3b.The inset shows a comparison with experimental data from the Constantine and Cain handbook 4 for large wt % and for T = 273 K (solid black line), 298 K (solid cyan line), and 323 K (black dashed line).Note that simulations correspond to T = 273 K (black circles), 293 K (cyan diamonds), and 323 K (black triangles-down).Fluctuations in the NPT ensemble 53 allow to compute the isothermal compressibility, κ T , the heat capacity, C P , and the heat capacity difference, C P − C V = TVα P 2 /κ T .Hence, it is possible to calculate the adiabatic compressibility from κ S = κ T C V /C P and compare it with experimental data.Note that C V is always lower than C P , otherwise it would imply κ T < 0 and a violation of the second law. 67For pure H 2 O 2 and T = 283 K, C V /C P yields 0.980 and decreases slightly with wt % at constant T. The model results for κ S are given in the main panel c of Figure 3, whereas the comparison with experiments 4 at T = 283 K and large wt % is given in the corresponding inset.The matching of both series is very good.
The dipoles gromacs tool computes the total dipole moment, M, of the simulation box as a function of time, its mean squared value, ⟨M 2 ⟩, and its mean value to the square, ⟨M⟩ 2 .From these quantities, it follows the relative dielectric constant, ϵ, 53,68 a property that impacts the penetration of the electric field into the liquid bulk and, hence, the long-range electrostatic interactions.The main panel d of Figure 3 shows the simulation results and the inset a comparison with , and the complete system.n h is given per number of oxygen atoms of H 2 O 2 molecules except for H 2 O -H 2 O, which is given in terms of the number of H 2 O molecules, and the total, given in terms of the total number of oxygen atoms.Insets without labeled axes share the same axes as the main panel.Experimental data in insets of (a) and (e) are taken from Easton et al. 62 and Giguere et al., 63 respectively.−66 Simulation outcomes are given in Tables S6−S13 of the SI.experimental data 4,66 for pure H 2 O 2 .Note that, despite performing relatively large runs, the data are noisy.We estimate error bars close to 5% (not shown to gain clarity).There is a general decreasing trend of ϵ with increasing temperature for all mixtures.However, the ϵ trend with wt % changes from increasing at low temperatures to decreasing at higher ones.In fact, for T = 313 and 323 K, the ϵ(wt %) curves are practically flat.Consequently, the ϵ window for the SPC/E model is smaller than for the H 2 O 2 one.We have found a very good agreement between simulations and the experimental results 4,66 for pure H 2 O 2 at 273 K, which worsens a little with increasing T. The agreement between our data and previous SPC/E results is also good. 69he enthalpy of vapor formation, ΔH vap , measures the overall attractive interaction between molecules, and it is frequently used to calibrate molecular models.It is given by where U g (T) is the potential energy of an isolated (gas) molecule at T, U l (T) is the potential energy of the bulk liquid per molecule at T, and R = 8.31446 J mol −1 .Note that this expression assumes that the specific volume of the liquid is much smaller than the corresponding gas and that the kinetic energy of vapor and liquid are the same.In addition, it does not account for quantum corrections.We can also assume that U g (T) = U g min + RT(3N sites − 6 − N c )/2, with N sites = 4, N c = 0, and U g min = −1.53kJ mol −1 for our H 2 O 2 model (the minimum potential energy), is the gas potential energy (the center of mass displacement and the rotation of the molecule around its principal axes do not contribute to the potential energy) 70 or compute U g (T) directly from an isolated molecule.In this last case, we constrain the displacement and rotation of the molecule so that only the degrees of freedom that interact with the internal potential energy are allowed.Both methods yield similar results. 71Note also that U g (T) is a constant for a rigid molecule such as the SPC/E.Our results compare well with the 48.8 kJ mol −1 obtained for the SPC/E model at 301 K. 72 Including self-polarization energy corrections 45 to the ΔH vap of SPC/E leads to 43.6 kJ mol −1 at the same T. 72 The ΔH vap results are shown in panel e of Figure 3.As expected, ΔH vap increases with wt % and decreases with T. In addition, the ΔH vap change with wt % increases with wt %.Our H 2 O 2 model and its mixtures with the SPC/E one follow the general trend of the experimental determinations of Giguere et al. 63 but overestimate ΔH vap for all wt % at T = 273 K.Note also that, as for SPC/E water, this overestimation may be explained in terms of the self-polarization correction.The experimental dipole moment of the isolated H 2 O 2 molecule is reported to range between 2.05 and 2.26 D (Constantine and Cain handbook; 4 Gross and Taylor 66 ).Considering the H 2 O 2 polarizability is 30% larger than that of water 73 and the 2.69 D value obtained for our model at 273 K, the self-polarization correction to the energy predicts a range of 6.6 to 3.0 kJ mol −1 .This range covers the difference between the predicted and experimental 63 ΔH vap values for pure H 2 O 2 .Finally, it is worth mentioning that, although we did not consider the ΔH vap as a target property, we took into account the vapor pressure and the surface tension, which are also strongly related to the intermolecular attraction.
As mentioned in the Simulations Details section, viscosities are calculated from nonequilibrium simulations 54 with an acceleration amplitude of 0.05 nm ps −2 .Results are shown in panel g of Figure 3.The μ values we have obtained for the pure SPC/E at T = 293 and 303 K are (0.756 ± 0.006) cP and (0.631 ± 0.005) cP, respectively, which can be compared with the value at 298 K of 0.729 cP obtained by Gonzaĺez et al. 77 (the experimental value is 0.896 cP 78 ).μ increases with H 2 O 2 concentration and decreases with T. The inset compares the simulation results with the experimental outcomes reported by Phibbs and Giguere 55 for pure H 2 O 2 .The disagreement with experiments diminishes with increasing T, but for T = 273 K, the overestimation reaches around 100%.In addition, viscosities were obtained by constraining the distance d OH and the angle a HOO .Without these constraints, the overestimation gets even worse.Thus, we can only claim that the μ trends are captured with our model.Note that the μ values given in Tables S1 and S2 frequently overestimate the experimental values. 55he last panel of Figure 3 shows the average number of hydrogen bonds per frame, n h , formed among H 2 O 2 molecules (leftmost inset), H 2 O molecules (left middle inset), H 2 O 2 with H 2 O molecules (right middle inset), and the total (rightmost inset).n h is given in terms of the number of oxygen atoms of H 2 O 2 molecules except for hydrogen bonds formed among H 2 O molecules, given per number of H 2 O molecules, and for the total, given in terms of the total number of oxygen atoms.The hydrogen bonds are computed by setting a cutoff distance for the donor−acceptor of 0.35 nm and a cutoff angle for the hydrogen-donor−acceptor sites of 30°.For the pure substances, n h is larger for H 2 O (1.86 at T = 273 K) than for H 2 O 2 (1.35 at T = 273 K).In addition, the n h ratio between the two quantities is practically constant when temperature increases from 273 to 343 K. On the other hand, n h among H 2 O 2 and H 2 O molecules reaches 2.24 at T = 273 K at high H 2 O 2 dilution.Overall, the total n h smoothly decreases with increasing wt %.Note that the n h dependence with T is small and similar for all hydrogen bonds.Note also that the increase of ΔH vap , λ, and T c with augmenting wt % is not due to the formation of more hydrogen bonds per oxygen atom but per molecule (from 1.86 to 2.7 at T = 273 K).Indeed, also the density of hydrogen bonds increases with wt %.In addition, the larger oxygen density also yields a larger LJ contribution.For T = 273 K, the ratio between n h for the pure substances is Journal of Chemical Theory and Computation 1.375 in favor of H 2 O, but the density ratio is 1.442 in favor of H 2 O 2 .Thus, n h in terms of volume increases 1.05 times when going from pure H 2 O toward pure H 2 O 2 , which partially explains the 1.16 ratio of ΔH vap .The former quantity slightly diminishes with increasing T, as also happens with the ΔH vap ratio.Finally, these results generally agree with the ones reported by Orabi and English, 18 although they report a larger n h for the pure H 2 O 2 substance (around 1.7 at 273 K).
Water−H 2 O 2 Liquid−Vapor Mixtures.As explained in the Simulations Details section, we access the vapor−liquid coexistence by placing a slab of liquid surrounded by its vapor in a prismatic cell with a height (z-axis) 3 times larger than its other two sides. 56,57This implementation is an alternative to the Gibbs ensemble technique, with the advantage that it allows accessing the surface tension (we do it following the virial route, 53 but one can also introduce area changes for this purpose 81 ).The obtained density profiles for H 2 O 2 , H 2 O, and the total are shown in panel a of Figure 4.The panel makes clear that these water and peroxide models produce fully mixed mixtures for all wt %.That is, there are no H 2 O 2 or H 2 O aggregates leading to bumps in the liquid profiles.There is, however, a slight preference for water to place at the liquid− vapor interface, displacing peroxide toward the bulk liquid.Indeed, the proportion of peroxide in the vapor phase is lower than in the liquid.
Panel b of Figure 4 is built from liquid−vapor profiles such as those shown in panel a.Here, the inner liquid−vapor coexistence corresponds to the pure SPC/E model and the outermost to the pure H 2 O 2 one.The critical points are obtained by using the law of rectilinear diameters 82 and a scaling law with a critical exponent of 0.32. 83,84For the SPC/E model, we get a critical point with T c = 622.8K, ρ c = 291.9kg m −3 , and P c = 126.6 bar, which is in good agreement with the ones previously reported. 85,86On the other hand, the critical point for the pure peroxide model leads to T c = (733 ± 8) K, P c = (224 ± 6) bar, and ρc = (446 ± 10) kg m −3 , which matches the experimental data reported by Nikitin et al., 87 T c = 728 K and P c = 220 bar (see the references therein for values obtained by other authors).In between the pure water and peroxide, the critical density and temperature smoothly vary as shown in Figure 4b.The inset of this same panel depicts the critical pressure against its temperature for the different peroxide concentrations.
For a fixed overall wt % of a simulation cell, H 2 O 2 (H 2 O) increases (decreases) its concentration in the bulk liquid and decreases (augments) it in the vapor phase.Thus, when plotting the vapor weight percentage of H 2 O 2 , wt % v , against the liquid one, wt % l , the points shift from the bisector of the axes to the right and bottom (see Figure 4c).The smaller the temperature, the larger the deviation from the wt % v = wt % l line.This behavior is in qualitative agreement with the experiments reported by Scatchard et al. 79 In addition, the wt % v (wt % l ) curves from experiments 4,79 and simulations practically coincide at low temperatures (see the black line and the circles obtained at T = 313 K of the inset).However, the good quantitative agreement worsens when increasing temperature (see the cyan line and the squares obtained at T = 393 K of the inset).Panel d of Figure 4 shows the surface tension, λ, as a function of wt %, whereas its inset provides a comparison with experimental data from Phibbs and Giguere 55 for T = 293 K.The surface tension values for the pure SPC/E model agree with previously reported data. 57,81The surface tension increases with wt % and decreases with T as a consequence of the varying strength of the intermolecular interactions, following a trend similar to that of ΔH vap .Thus, it is not very surprising that our H 2 O 2 model overestimates the λ value.At T = 293 K, our value is (86.3 ± 1.2) dyn cm −1 , which is 7% larger than the 80.5 dyn cm −1 value reported by Phibbs and Giguere. 55However, the overestimation disappears for 80 wt % and changes sign for wt % < 80, as the SPC/E model underestimates the water surface tension.
Finally, the vapor pressure is depicted in panel e of Figure 4 as ln(P v ) as a function of 1/T.This representation is common due to the Clausius−Clapeyron equation, ln(P v1 /P v2 ) = ΔH vap (1/T 2 − 1/T 1 )/R, from where ΔH vap can be computed.However, the error in the determination of P v does not allow for a good estimation of ΔH vap through this path.Indeed, the relative error of P v strongly increases with decreasing T when employing the virial route.Thus, for P v < 2 bar, we are computing P v from ρ v by assuming an ideal behavior for the vapor phase.In turn, these data are employed to estimate the boiling temperature by fitting a second-order polynomial to P v (T) for low P v and its intersection with P v = 1.013 bar.The resulting boiling temperature is depicted against wt % and compared with experimental data from Giguere and Maass 80 in the inset.Note that the agreement is good when close to 100 wt %.

■ CONCLUSIONS
−28 The model does not include an LJ contribution to hydrogen sites and includes an RB dihedral function.The ANN optimization allowed us to consider several thermodynamic properties of water−hydrogen peroxide mixtures evaluated at different temperatures and compositions as targets.We explored the results of several ANNs and selected the ones that minimized the differences between computations and experimental data.For this last step, we assigned weights to each property based on the uncertainties of the simulation results and composition.We set lower weights when considering mixtures to avoid a strong dependence on the selected water model (SPC/E 45 ).However, we believe it is important to target not only properties for pure substances but also mixtures.
2][63][64][65][66]79,80 The best matches are found for properties such as density as a function of temperature, adiabatic compressibility, relative dielectric constant, diffusion coefficient for diluted mixtures in water, critical properties, and vapor composition at low temperatures. Howevr, there are larger mismatches observed for viscosities at low temperatures and vapor composition at high temperatures, which are overestimated by the model.We believe that the model predictions for mixtures could improve with the consideration of more sophisticated water models such as TIP4P, 24 TIP4P-2005, 88 or TIP4P-ϵ.25 Note that creating a model that can fit all properties under all conditions is a challenging task.For example, decreasing the attraction of the LJ parameter between oxygen atoms could improve the enthalpy of vapor formation but worsen the vapor composition.Furthermore, the interplay between all parameters is not trivial.Fortunately, ANNs can take this complex interplay into account and may become the preferred tool for designing new molecular models.
The itp file defining the H 2 O 2 model; tables to train (Table S1) and test (Table S2) the artificial neural networks (ANNs), results from different ANNs (Table S3), optimized parameters from the ANNs (Table S4), and the final results from a local search around the outcomes of the ANNs (Table S5).Tables S6

Figure 1 .
Figure 1.(Top left) Scheme showing the geometrical parameters of the H 2 O 2 model.Side (bottom left), top (bottom right), and front (top right) views of the ball-and-stick representation of the model.The labels d OO and d OH correspond to the oxygen−oxygen and oxygen−hydrogen distances, respectively.The a HOO and a d labels refer to the hydrogen−oxygen−oxygen and dihedral angles, respectively.
constraints c 0 = 0.312c 2 , c 1 = −c 2 , c 3 = 0.300c 2 , and c 4 = c 5 = 0.With these constraints and for c 2 ranging from 25.5 kJ mol −1 to 36.0 kJ mol −1 , we get a function with minima at a d ≈ ±130°, a large maximum at a d = 0°, and a tiny barrier at a d = ± 180°(see the cyan line inside panel b of Figure

Figure 2 .
Figure 2. (a) Radial distribution function, g(r), for the O−O (solid black line), O−H (black dashed line), and H−H (solid cyan line) pairs (note that the peaks corresponding to the internal O−O and O−H bonds are not considered to gain clarity).The inset shows the a HOO probability distribution function, pdf.(b) Probability distribution function, pdf, for the dihedral angle a d (black) and the Ryckaert-Bellemans potential contribution (cyan line and right y-axis).The inset depicts the probability distribution function of H 2 O 2 dipole moments.All data correspond to a bulk liquid of pure H 2 O 2 , T = 293 K, and P = 1 bar.
Finally, we set the single bond O−O force constant equal to that of the N−O single bond of the OPLS-aa force field and the H−O force constant equal to the one of the TIP4F 46 water model.It should be noted that there is no O−O single bond defined in the OPLS-aa force field, but the C−O and N−O bonds have the same force constant values due to the similarity in size and electronegativity between C and N.As O is also similar to N and C in size and electronegativity, we expect the O−O bond to have a similar strength.

Figure 3 .
Figure 3. (a) Density, (b) thermal expansion coefficient, (c) adiabatic compressibility, (d) relative dielectric constant, (e) enthalpy of vaporization, (f) diffusion coefficient, (g) viscosity, and (h) the number of hydrogen bonds per oxygen atom as a function of the weight percentage of H 2 O 2 in the water−peroxide mixture.For all curves, circles, squares, diamonds, triangles (up, left, down, right), and plus symbols correspond to T = 273, 283, 293, 303, 313, 323, 333, and 343 K, respectively.The insets compare the simulation results (symbols) with experimental data (lines) for (a) T = 273 and 313 K, (b) 273, 298, and 323 K, (c) 283 K, (d) pure H 2 O 2 , (e) 273 K, and (f) pure H 2 O 2 .The main panel (f) shows the diffusion coefficient for H 2 O 2 molecules and the inset for H 2 O.The four plots of panel (h) depict, from left to right,n h for H 2 O 2 -H 2 O 2 , H 2 O -H 2 O, H 2 O 2 -H 2 O, and the complete system.n h is given per number of oxygen atoms of H 2 O 2 molecules except for H 2 O -H 2 O, which is given in terms of the number of H 2 O molecules, and the total, given in terms of the total number of oxygen atoms.Insets without labeled axes share the same axes as the main panel.Experimental data in insets of (a) and (e) are taken from Easton et al.62 and Giguere et al.,63 respectively.Experimental data for all other insets are taken from the Constantine and Cain handbook,4 which summarizes and provides fits from data of different sources.55,64−66Simulation outcomes are given in Tables S6−S13 of the SI.

Figure 4 .
Figure 4. (a) Density profiles for H 2 O 2 (solid black lines), H 2 O (solid cyan lines), and total (black dashed lines).All curves correspond to T = 533 K, and the wt % increases from left to right and from top to bottom (10, 20, 30, 40, 60, and 80).(b) Vapor−liquid phase diagram.Critical points are shown as unconnected symbols.The inset shows the critical pressure against the critical temperature.(c) wt % v against wt % l of the coexisting phases.The inset compares simulation results for T = 313 K (circles) and 393 K (squares) with experimental data 79 (black lines for T = 313 K and cyan lines for T = 394 K).(d) Surface tension as a function of wt %.The inset compares the outcomes from the model (circles) with experimental data 55 (solid line) for T = 293 K. (e) Natural logarithm of the vapor pressure against the reciprocal temperature.The inset compares the predicted boiling point temperatures (circles) with the experimental values 80 (solid line) as a function of wt %.In panels (b) and (e), circles, squares, diamonds, triangles (up, left, down, right), and plus symbols correspond to 0, 10, 20, 30, 40, 60, 80, and 100 of wt %.In panels (c) and (d), circles, squares, diamonds, triangles (up, left, down, right), and plus symbols correspond to T values from 293 to 433 K in steps of 10 K. Insets without labeled axes share the same axes as the main panel.Simulation outcomes are given in Tables S14 and S15 of the SI.